function VARIAN = HPZ_Varian_efficiency_index (expenditure)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
% implementing the algorithm in Varian (1993) using the version in Alcantud
% et al (2010) (Algorithm 3). Note that the example in this paper is wrong.

global Threshold

[rows,cols] = size (expenditure);

var = ones(rows,1);

GARPE = 1;

while GARPE
    
    % Line 4
    exp_var = expenditure - diag((diag(expenditure)).*(ones(rows,1)-var));
    
    %The matrix REF has at the cell in the i'th row and the j'th
    %column, the difference between the value of the bundle that was chosen in 
    %observation i and the bundle that was chosen in observation j given the 
    %prices of observation i

    REF = diag(exp_var)*ones(rows,1)' - exp_var;
    
    %The matrix RATIO has at the cell in the i'th row and the j'th
    %column, the inverse of the ratio between the value of the bundle that was chosen in 
    %observation i and the bundle that was chosen in observation j given the 
    %prices of observation i
    
    RATIO = exp_var./(diag(exp_var)*ones(rows,1)');
    
    %The matrix DRP has at the cell in the i'th row and the j'th
    %column, 1 if and only if the bundle that was chosen in 
    %observation i is directly reveal prefered to the bundle that was chosen 
    %in observation j (the corresponding value in REF is greater than -1).

    DRP = ceil((REF+Threshold)/(max(max(abs(REF+Threshold)))+1));
    
    %The matrix SDRP has at the cell in the i'th row and the j'th
    %column, 1 if and only if the bundle that was chosen in 
    %observation i is strictly directly reveal prefered to the bundle that was chosen 
    %in observation j (the corresponding value in REF is greater than -1).

    SDRP = ceil((REF-Threshold)/(max(max(abs(REF-Threshold)))+1));

    % statement needed for the graph theory external package

    set_matlab_bgl_default(struct('full2sparse',1));

    %The matrix NS_RP has at the cell in the i'th row and the j'th
    %column, Inf if and only if the bundle that was chosen in 
    %observation i is not reveal prefered to the bundle that was chosen 
    %in observation j. Otherwise it includes a positive integer.

    NS_RP = all_shortest_paths(DRP);

    %Create RP

    RP = eye(rows,cols);

    %The matrix RP has at the cell in the i'th row and the j'th
    %column, 1 if and only if the bundle that was chosen in 
    %observation i is reveal prefered to the bundle that was chosen 
    %in observation j. 

    for j=1:rows
        for k=1:cols
            if ~isinf(NS_RP(j,k)) && ~(k==j)
                RP(j,k)=1;
            end
        end
    end

    %To test for GARP we will do the following: for every pair of
    %choices x and y if xRy then not yP0x. We will take RP and the transpose of
    %SDRP and multiply element by element. Every 1 correspondes to 
    %xRy and yP0x. The final matrix is the zero matrix if and only if 
    %GARP is satisfied. 

    GARP = RP.*(SDRP');

    GARP_ERRORS = sum(sum(GARP));

    %If the GARP matrix is only zeros then the data satisfies GARP
    if GARP_ERRORS==0
        GARPE=0;        
    else
        % Line 7 - Gv(x_j) is the jth column of GARP
        
        % Line 8 - Pert is a row vector. if the jth column is zero then Gv(x_j) is empty. 
        % otherwise it includes the number required in line 8.
        Pert_mat = GARP.*(RATIO');
        
        for j=1:rows
            for k=1:cols
                if Pert_mat(j,k)>=1
                    Pert_mat(j,k)=0;
                end
            end
        end
        
        % Line 9
        r = max(Pert_mat);
        
        [v,w] = max(r);
        
        % Line 10
        var(w) = v * var(w);
    end    
end

VARIAN = [min(var),mean(var),sum(var.^2)];

end

